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We obtain thermostatted ring polymer molecular dynamics (TRPMD) from exact quantum dynamics via 
Matsubara dynamics, a recently-derived form of linearization which conserves the quantum Boltzmann 
distribution. Performing a contour integral in the complex quantum Boltzmann distribution of Matsubara 
dynamics, replacement of the imaginary Liouvillian which results with a Fokker-Planck term gives TRPMD. 
We thereby provide error terms between TRPMD and quantum dynamics and predict the systems in which 
they are likely to be small. Using a harmonic analysis we show that careful addition of friction causes 
the correct oscillation frequency of the higher ring-polymer normal modes in a harmonic well, which we 
illustrate with calculation of the position-squared autocorrelation function. However, no physical friction 
parameter will produce the correct fluctuation dynamics for a parabolic barrier. The results in this paper 
are consistent with previous numerical studies and advise the use of TRPMD for the computation of spectra. 
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I. INTRODUCTION 

The computation of thermal time-correl atio n functions 
is of central importance in chemical physicJ^^^ in order to 
evaluate many physically observable quantities such as 
reaction rates, diffusion constants, spectra and scattering 
data.l^ 

Exact evaluation of the quantum correlation function 
scales exponentially with swtem size and so is impractical 
for more than a few atoms.^ There is consequently a need 
for computationally tract able approximations to quantum 
time-correlation functiond^, preferably which are known 
to be equivalent to the quantum result in certain limits, 
and for which the likely error is known in advance of 
calculation. One crude solution is to use a purely classical 
correlation function, which will scale linearly with system 
size. However, the classical Boltzmann distribution is 
highly inaccurate for many systems such as water at 
room temperatur^, a nd i gnores effects such as tunnelling 
and zero-point energy.I^ There is consequently a need to 
incorporate quantum statistics into such calculations, but 
with approximate, preferably classical-like, dynamics. 

Various approaches have been developed, including the 
“classical Wigner” or linearized semiclassical initial value 
representation (LSC-IVR) methocP^^HH which truncates 
the (exact) Moyal series^ for time evolution at h^, Cen¬ 
troid Molecular Dynamics (CMDjPiHi^, which propagates 
the path-integral centroid in the mean-field of the other 
path-integral normal m odes, a nd Ring Polymer Molec¬ 
ular Dynamics (RPMDjPE^^^ which takes the classical 
dynamics of a ring polymeil^ literally. 

All these methods have various limitations; LSC-IVR 
does not conserve the quantum Boltzmann distribution 
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leading to zero-point energy leakag^®, where as CM P and 
RPMD both fail for multidimensional spectr£p21241 CMD 
has the curvature problem where peaks are broadened 
and red-shifted whereas RPMD has spurious resonances 
where the ring polymer springs couple to frequencies in 
the potential leading to splitting of the physical peak.^^ 
Recently, Thermostatted Ring Polymer Molecular Dy¬ 
namics (TRPMD) has been introduced, which applied 
a Langeyint^rmostafpSH^ to the dynamics of the ring 
polymeJ ^^ I ^^ I ^° l This was origmally conceived for the 
evalutation of static properties^, but it appeared to be 
remarkably successful for the computation of spectrsP^*^, 
accurately replicating multidimensional spectra where 
CMD and RPMD fail and correctly predicting the dif¬ 
fusion and rotational constants of liquid water. Like 
RPMD, the short-time, transition-state theory (TST) 
limit of the TRPMD flux-side correlation functio n is id en- 
tical to quantum transition-state theory (QTST^pSEHMI: 
the instantaneous thermal ring-polymer flux through a 
position-space dividing surface is equal to the intanta- 
neous thermal quantum flux, and the TRPMD rate will 
equal the exact quantum rate in the absence of recrossing 
by either the quantum dynamics or TRPMD dynamicp21_ 
Because TRPMD obeys detailed balance, its reaction rate 
is independent of the location of the dividing surfac^^, as 
is the case for RPMD and CMD but not many TST-based 
methods. 


Nevertheless, TRPMD is not without its faults; like 
RPMD and CMD it fails to capture effects such as a 
Fermi resonance involving a fourth-order coupling in the 
Zundel catiorP^, and beneath the crossover temperature 
(see Eq. (29)) application of friction to reaction rates 
causes them to decrease, resulting in a less accurate result 
compared to RPMD for symmetric systems, and a more 
accurate result (but with an adjustable parameter whose 
value is not determinable in advance) for asymmetric 
systems.*^ 
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Very recently, both RPMD and CMD have been ob¬ 
tained from the exact quantum time-correlation func¬ 
tion via “Matsubara dynamics”, a form of lineariz ation 
which conserves the quantum Boltzmann distributioiP^^^. 
Matsubara dynamics results from discarding fluctuations 
of the very high frequency path-integral normal modes 
(higher frequencies than those required to converge the 
quantum Boltzmann distribution) from the (exact) Moyal 
bracket and is inherently classical as well as satisfying 
detailed balanc^^. However, Matsubara dynamics is not 
amenable to computation in large systems since it suf¬ 
fers from the sign problem due to a phase factor in the 
complex quantum Boltzmann distributiorP^. An approxi¬ 
mation to Matsubara dynamics where the centroid moves 
in the mean field of the other Matsubara modes leads to 
CMD, and by moving the momentum contour in the quan¬ 
tum Boltzmann distribution and discar ding the imaginary 
Liouvillian which results, RPMD ariseP^. 

Obtaining RPMD and CMD from the exact quantum 
expression provides analytical expressions for their error 
from the quantum result, such that it can be known a 
priori whether they will function well in a given system, 
whereas previously one had to rely on induction from 
earlier numerical studies on systems for which RPMD or 
CMD had been successful, though there was no guaran¬ 
tee that such reasoning would extend to a new system. 
This was seen, for example, in RPMD rate theory failing 
in the Marcus inverted regim^^ despite being very suc¬ 
cessful for rate computation in a large variety of other 
gygtempSmiHlT] derivation of QTSlI^SIlll, it 

can be known a priori that a system whose optimal ring- 
polymer dividing surface will be significantly recrossed by 
the quantum dynamics will not have its rate accurately 
computed by RPMEP^. 

Consequently, investigating whether TRPMD could 
also be obtained from exact quantum time evolution and 
thereby discerning the situations where it is likely to work 
a priori, rather than relying on the (small but g rowing) 
literature of its application to physical systemP^EoM]^ 
would be of considerable benefit to the field. 

In this paper we obtain TRPMD from exact quantum 
dynamics by showing that it is a stochastic approximation 
to Matsubara dynamics. To obtain a computationally 
tractable approximation to Matsubara dynamics, we move 
the momentum contour in the complex plane in order 
to convert the complex quantum Boltzmann distribution 
into the real ring polymer distribution. This transforma¬ 
tion generates a complex Liouvillian in the dynamicJ^, 
which is not in itself amenable to computation due to 
the complex trajectories which result. Previously, the 
imaginary part of the Liouvillian was simply discarded, 
shifting the oscillation frequencies in the higher normal 
modes and leading to RPMD, but here we replace it with 
a Fokker-Planck term, producing TRPMD. 

To examine the effect of the friction matrix we conisder 
a harmonic well and a parabolic barrier, for which the cor¬ 
relation function (if defined) can be evaluated analytically. 
We find that a unique and system-independent value of 


the friction matrix causes all normal modes to oscillate at 
the correct (external) frequency in a harmonic well, and 
illustrate this with the position-squared autocorrelation 
function, where TRPMD has the correct zero-time value 
and frequency; neith er RP MD nor CMD can reproduce 
both these propertie^^^EH then examine a parabolic 
barrier, where CMD and RPMD have the incorrect fluc¬ 
tuation dynamics; the higher normal modes in RPMD 
being bound (above the relevant crossover temperature, 
see Eq. rather than unbound as in Matsubara dy¬ 

namics. Here application of any meaningful (i.e. positive) 
friction does not cause the erroneously bound normal 
modes in TRPMD to become scattering, and nor does 
it cause unbound modes to have the correct escape fre¬ 
quency, meaning that application of friction is unlikely 
to assist in the accuracy of reaction rate or diffusion 
calculation.!^ 

We begin by revisiting Matsubara dynamics in sec¬ 
tion]^ followed by obtaining TRPMD in section m and 
examining the friction matrix in section IV, before pre¬ 


senting conclusions in section VI 


II. SUMMARY OF MATSUBARA DYNAMICS 


For simplicity, we consider a one-dimensional system 
with mass m, co-ordinate q, and Hamiltonian H = 
p^/2m + V{q), where V{q) is the potential.!^ Extensions 
to further dimensions follows immediately and merely 
requires more indices. The Kubo-transformed thermal 
quantum time-correlation function at inverse temperature 
P = l/k^T ipl 

/ X 1 

CAB{t) = - 

( 1 ) 


daTr 




and which can easily be related to the conv entio nal 
asymmetric-split quantum correlation functioiPE^. If 
A and B are linear operators in position or momen¬ 
tum, Eq . (|T1) is identical to the Generalized Kubo 

transfornPlte^ 


N-1 

X (gi_i - Ai_i/2|e“^*|gi-h Ai/2) 

1=0 

X + - A,/2) 

( 2 ) 

where / dq = ^ dqi (similarly for f dz and 

JdA), 


^(q) 


1 

N 


N-l 

E ^(9*) 


2=0 


(3) 


and likewise for B{q). 
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The full derivation of Matsubara dynamics is given in 
Ref. [36! and here we sketch the relevant details. To obtain 
the time-evolution of Eq. (H in the phase-space represen¬ 
tation, we take its Wigner Transf orirP and differentiate 
w.r.t. time to obtain a Moyal serieJ^^lSIl which is formally 
exact. We then transform the correlation function from 
bead co-ordinates (p, q) to ring-polymer normal mode 
CO- ordinate^ (P, Q), as detailed in appendix [A] such 
that Qo and Pq are the position and momentum centroids 
respectively. 

Truncating the resulting Moyal series (either in the 
normal mode or bead representation) at 0{hP) leads to 
the lineari zed sem iclassical initial value representation 
(LSC-IVR)pl2li2EHl^ which involves propagating trajecto¬ 
ries under the classical Hamiltonian of the system drawn 
from a Wigner-transformed quantum Boltzmann distri¬ 
bution. Conversely, truncating the Moyal bracket to the 
lowest M ‘Matsubara’ normal modes (see appendix 
results in a dynamics which is inherently classical (all 
powers of and higher vanish from the Moyal series 

without further approximation) and which conserves the 
quantum Boltzmann distr ibut ion and satisfies detailed 
balancd^, unlike LSC-IVRP^. This leads to a classical- 
like correlation functiorP^ 




= ^ [ dP [ dQ e-/5[^M(P,Q)-*eM(P.Q)] 

2Trh J J 

X H(Q)e^Mi‘s(Q) (4) 


computation in complex systems. To make the distri¬ 
bution real, we continue into the complex plane of P 
with 


Pj = Pj — imuijQ-j 


( 8 ) 


for all j (such that no analytic continuation is necessary 
for the momentum centroid) to give 


r[^^] _ r[^] j_ ^ 


(9) 


where is the ring polymer Liouvillian, 


ri^] _ 

^RP “ 


(M-l)/2 

E 

3 = -{M-\)j1 


PiJL 

m dQj 


ai7[^i(Q) 

dQj 


mujjQj 


and the imaginary component of the Liouvillian is 

(M-l)/2 


j=-iM-l)/2 


~ , r. d ^ d 

1 Pn ^ ^ ^ j ' 


'dP-, 


dQ-, 


dPj 

( 10 ) 


( 11 ) 


This transformation also converts the complex Matsubara 
distribution into the real ring polymer distribution. 


g-/3[ffM(P,Q)-*e„(P,Q)] ^ g-/3flM(P,Q) 


( 12 ) 


where the Matsubara Hamiltonian is 


where the ring-polymer Hamiltonian is 


(M-l)/2 p2 

HMiP,Q)= E ^ + 

t7['^l(Q) is defined in the appendix, um = — 

l)/2]!^, and the phase factor is 

(M-l)/2 

0m(P,Q)= 5] (6) 

i=-(M-l)/2 


where Coj = 2TTj/j3h are the Matsubara frequencies 
which, in this definition, can be negative. The integrals 
are taken to mean f dP = ^Pj ^^*1 

likewise for J dQ. Matsubara dynamics is defined by the 
Liouvillian 


such that = {•, L7m(P, Q)} where {•, •} is the clas¬ 
sical Poisson bracket .1^ 


III. EMERGENCE OF TRPMD 

The Matsubara correlation function in Eq. suffers 
from the sign problem, such that it is not amenable to 


Rm(P, Q) 


(M-l)/2 

E 

i=-(M-l)/2 



+[/[“! (Q). 
(13) 


Both and independently conserve the quantum 
Boltzmann distribution. 

In AppendixO we prove that the complex dynamics 
generated by Cp is analytic everywhere in the complex 
plane, and by contour integration of Eq. Q it rigorously 
follows 


^[M], . _ aM 


dP J dQ 

X H(Q)e''^kp' 


(4*?+*4"')*i3(Q) 


Sit) 


(14) 


where S (t) corresponds to the vertical edges of the inte¬ 
gration contour; in Appendix we give evidence to show 
that in many cases the edge term will vanish, though 
for an arbitrary system propagated to a finite time it is, 
strictly speaking, part of the error term between Matsub¬ 
ara dynamics and RPMD/TRPMD. 


Although the real ring-polymer distribution in Eq. (14) 


would, prima facie, allow evaluation of the correlation 
function by inexpensive Monte Carlo techniques, the pres¬ 
ence of in Eq. ([I4|) causes unstable trajectories to 
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which are no easier to treat numerically than 
the sign problem in the complex Matsubara distribution!^ 
In previous researclP^ it was shown that approximating 
Eq. (141 by discarding in order to make the tra¬ 


jectories real but still conserve the quantum Boltzmann 
distribution, produces RPMD. This approximation raises 
the oscillation frequency of the higher (j ^ 0) normal 
modes; in a harmonic potential with external frequency 

ujh they oscillate at uij = This is the origin of 

the ‘sp uriou s resonances’ of RPMD in multidimensional 
spectrd;23124] qualitative failure of RPMD at calcu- 

lating the position-squared autocorrelation functiorP^IlIl 
It was also shown that a mean-field approximation to 
Eq. (14), where the centroid moves in the mean field of 
the other ring polymer modes, leads to CMEI^. 

This naturally motivates investigating whether there 
is some other approximation to the dynamics in Eq. (14) 


which, like RPMD, is real and conserves the quantum 
Boltzmann distribution, but unlike RPMD has the correct 
oscillation frequencies of the higher normal mode^^, and 
possibly also has the correct fluctuation dynamics at 
barriers. Addition of a friction (Langevin) term to the 
dynamics of a harmonic oscillator is known to decrease the 
oscillation frequencjES and we therefore define a stochastic 
dynamics by the Eokker-Planck adjoint operator 


,[M]t _AM\ -[M]t 

■^RP “•‘•'RP ' -^wn 


RP 


(15) 


where the white-noise thermostat which conserves the 
ring-polymer distribution is 


= -p. r. Vp + -Vp • r. Vp ( 16 ) 
with r an M X M positive semidefinite friction matrix. 


This allows us to approximate Eq. (14) as 


^[M] 


<''-S 


dP / dQ e 




^U(Q) 


e-^Rp 


‘B(Q), 

(17) 


Matsubara dynamics, we seek to determine an optimal 
friction parameter to reduce or remove the ‘side-effects’ 
of the analytic continutation to form RPMD, namely the 
incorrect frequencies of the higher normal modes in a 
bound systerrP, and the incorrect fluctuation dynamics 
in an unbound (scattering) systenP^. 


A. Harmonic well 

Eor a model bound system, we study the harmonic 
potential 


V{q) = ^mujlq^ (18) 

for which the ring polymer normal modes decouple and 
the dynamics can be solved exactly, and we detemine 
which elements of a diagonal friction matrix will cause 
oscillation at a correct (external) frequency ujh, as is 
the case for analytically continued Matsubara dynamics 
[Eq. ( [l^ ] in a harmonic potential,!^ 

P LJ 

QAt) = Qn cos(uJht) H- — siniujht) + i — Q-, sin(a;;jt). 

muJh 

(19) 

For TRPMD, the trajectories are not deterministic and 
we define the time-evolved phase-space density Qj{t) = 
Qj{Qj,Pj,t) which is evolved with from initial 

conditions of {Qj, Pj) at t = 0. For a harmonic potential 
and moderate (underdamped) friction , the time-evolution 
of Qj{t) can be solved analytically as^ES 


Qj{t) = e 


Qj cos(wjt) 


+ 


) sin(wA) 

mujj 2u}j ' ^ ^ ’ 


( 20 ) 


where the observed (damped) frequency of oscillation is 


which is TRPMD.l^The error term between the quantum 
result and TRPMD is therefore discarding the dynamics 
of the highest {N — M) normal modes to give Matsubara 
dynamics (see Eq. (B2) of Ref. [36]), the edges of the 
contour used in analytic continuation (which we suspect 
to be zero, see Eq. (CIO)), and the difference between the 
TRPMD and Matsubara propagators, namely — 

j[M]t 

vHwn 


IV. FRICTION CONSIDERATIONS 

There are already numerical studies of the effe ct of 
friction on various quantities computed by TRPMEPHSH^ 
and here we take a more theoretical approach in light 
of obtaining TRPMD from quantum dynamics in the 
previous section. Since TRPMD is an approximation to 


V“2 + 'd-rt/'i- 


( 21 ) 


We immediately see Tjji = 2\ujj\Sjji will ensure that 
LOj = uJh and oscillation at the correct external frequency, 
a result previously suggested on the grounds of minimizing 
the Hamiltonian correlation time for a ring polymer in 
a harmo nic po tential, and thereby optimizing statistical 
sampling2S122l_ Xo investigate different friction strengths 
related to this we therefore define a parameter A such 

Although the position-squared and momentum-squared 
correlation functions will oscillate at the external fre¬ 
quency with A = 1 (see below), examination o f the 
position-position spectrum for a given normal mod^^sUZl 


^TRPMD 


(cu) oc 


1 

{ujI + - Cli2) -I- ry2^2 


( 22 ) 
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shows that the maximum in the spectrum will be at a; = 
\l^h + ‘^1 — r|^/2, suggesting a friction parameter of 

A = Furthermore, consideration of the momentum 

spectrum = m?{ oj) shows that 

the maximum in the momentum spectrum is always at the 
(erroneously high) ring polymer frequency ui = 
and increasing friction merely broadens the peak. 


B. Numerical example 


exactly equal to the quantum result but both RPMD and 
CMD fail to qualitatively reproduc^^. CMD produces 
the incorrect result at t = 0 but then oscillates at the 
correct frequency (though incorrect amplitude), whereas 
RPMD is exact at zero time but then deviates wildly from 
the quantum result at finite time due to the prese nce of 
the spurious frequencies in the higher normal mode^^^liS] 
The “nonlinear operator” problem for which this is the 
archetypal model not only occurs in toy s ystem s but is 
also observed in inelastic neutron scatterin^^Ml 

The exact quantum position-squared auto corre lation 
function in the harmonic potential Eq. (18) i^SUlll 


CMD, RPMD and TRPMD correlation functions and 
spectra h ave already been the subjec t of many numeri¬ 
cal whose results have 

been broadly summarized in the introduction. To clar¬ 
ify the nature of the approximations inherent in CMD, 
RPMD and TRPMD from Matsubara dynamics, we ex¬ 
amine the position-squared autocorrelation function for 
a harmonic oscillator, for which Matsubara dynamics is 


Cq2q-2 (t) — 








coth 




cos(2w/i<) 


2coth^ 1 


(23) 


which in appendix]^ we show is exactly replicate d by the 
Matsubara correlation function. For RPMD, it 


^RPMD 


w = 


1 


(M-l)/2 

E 


1 


2cos^[(a;^ + 




- w" 


(M-l)/2 

E 


' i=-(M-l)/2 ^ J ri j k^-{M-l)/2 

whereas the TRPMD result for the optimal damping frequencies Pjj/ = 2\ojj\6jj' is 




/-fXRPMD 


(M-l)/2 

w=E 


1 


2g-2|(i3|t 


- Co] 1 -h Co] 


UJj 


cos{iOht) -\ - - sm{iOht) 

OOh 


(M-l)/2 

E 

fc=-(M-l)/2 


-I- Cbj, 


(24) 


(25) 


For comparison, the CMD position-squared autocorre¬ 
lation function (using the CMD with classical operators 
inetho(PE32lI13EIl) ig 

C'p^2°(t) = [2 cos{iOhtf + 1] . (26) 

We use parameters to facilitate comparison with previous 
literatur^^ h = k^=m = uoh = ^ and results for 
systems of varying /3 are presented in Fig. 

At high temperatures (/? = 1), all methods are a good 
approximation to the quantum result and the RPMD and 
TRPMD results are indistinguishable to within graphi¬ 
cal accuracy. At /3 = 4, the amplitude of oscillations is 
incorrect for all methods, though TRPMD starts at the 
correct value whereas CMD is too low. The RPMD corre¬ 
lation function shows deviations from harmonic behaviour 
due to the higher normal modes. At /? = 10, the CMD 
correlation function is far too small and RPMD cannot 
replicate the oscillations, unlike TRPMD. 

We then examine the effect of different friction param¬ 
eters in Fig. choosing the /3 = 10 system to exemplify 


the effect of damping. The A = 0 (RPMD) result oscil¬ 
lates erratically, as in the third panel of Fig. E Applying 
very small friction (A = 0.1) noticeably improves the 
correlation function but contamination from higher nor¬ 
mal modes is still evident. Around the optimal damping 
(A = 0.5, A = 1 and A = 1.5) the correlation functions 
are extremely similar, settling to the correct frequency 
(though incorrect amplitude) after one oscillation. In¬ 
creasing the friction yet further (A = 5) causes the correct 
oscillation frequency (as all modes apart from the cen¬ 
troid are overdamped) but the slow decay of the heavily 
overdamped higher normal modes causes midpoint of the 
oscillation to decay slowly over time. These results (which 
can be derived analytically for the harmonic oscillator) 
are broadly consistent with those observed numerically 
in more complex systems (see Fig. 3 of Ref. [25]l . where a 
broad range of friction parameters around 0.5 < A < 1.5 
led to similar results. 

For the position-squared autocorrelation function at 
least, it seems that TRPMD with A = 1 combines the best 
features of both RPMD and CMD; the correct zero-time 
value and the correct amplitude of oscillation, and that 
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C. Parabolic barrier 


For an unbound, scattering system, or where barrier 
dynamics are required such as thermal rate calculation, 
we instead consider a parabolic barrier with a potential 

V{q) = -lmujl (27) 

We firstly observe that RPMD (and CMD) have qualita¬ 
tively incorrect fluctuation dynamics at barriers“ while 
all the Matsubara modes are scattering. 


- 

QJt) cosh(uJbt) H- — sinh(wbf) 

mwb 

+ (5_,-sinh(a;bt) (28) 

w 

the RPMD higher normal modes are generally bound, with 
a frequency of Wj = yjCoj — . As the temperature is 

lowered, modes become successively unbou nd, b eginning 
with j = ±1 at the ‘crossover’ temperatur^^UllI 


/3c = 


27r 

hiOb 


(29) 


and the jth normal mode will become unbound when 
P > \j\Pcy but with a scattering (imaginary) frequency 

of yjujb — Cjj. Despite these shortcomings, RPMD has 

been very accurate for rate calculation, partly because 
of the correct short-time TST limilP^ (see Introduction), 
and since reaction above the crossover temperature is 
dominated by the motion of the centroid, which scatters 
at the correct imaginary frequency in RPMD, TRPMD 
and CMD. 

Considering the jth mode at a temperature /3 < \j\Pc 
such that it is bound in RPMD {ojb < |wj|), very weak 
friction {Tjj/2 < uij — uj^) leads to damped oscillatory 


motion as in Eq. (201, but with u}j = yj— 

Stronger friction leads to an overdamping solution. 


Qj(t) = e 


Qj cosh(Cjt) 


FIG. 1. Position-squared autocorrelation function for a har¬ 
monic oscillator, with j3 = 1 (top), j3 = A (middle) and /3 = 10 
(bottom). Black circles, quantum; solid green line, TRPMD 
(with optimal damping); red dots, RPMD; blue dot-dashes, 
CMD. 


there exists a sizeable range of friction parameters around 
A = 1 in which these qualitative features are captured. 
However, TRPMD with optimal friction is by no means 
perfect; the amplitude of oscillations decays within one 
oscillation and is far too small, since by this time all 
modes apart from the centroid are essentially completely 
damped. 


-b 


mCj 


Qj^n 

2Q 


sinh(Cjt) 


(30) 


where Q = — uj'j + {Tjj/2Y can be considered the 

imaginary frequency counterpart to Wj. 

The presence of the prefactor in Eq. (301, which 


in the oscillatory case of Eq. (20) causes damping but 


leaves the frequency untouched, means that no physical 
(i.e. real and positiv^^ value of the friction parameter 
exists which would make Eq. ( |30[ ) have an unbound so¬ 
lution. Increasing friction merely causes the oscillator to 
become more overdamped. 

If the normal mode is unbound in RPMD [fi > \j\Pc 
and therefore ujb > |wjD then Eq. (30) still holds, but 
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Quantum 
A = 0 
A = 0.1 
A = 0.5 
A = 1 
A = 1.5 
A = 5 


FIG. 2. Position-squared autocorrelation function for a harmonic oscillator at /I = 10, showing the exact quantum result and 
TRPMD at a varying friction parameters. For clarity, the figure is zoomed in around the TRPMD correlation function. 


the solution has a scattering component for all Tjj since 
—Tjjl2 -h Qj > 0. We would like to increase the escape 

rate from the RPMD value of to Wh, but adding 

friction only decreases the rate of escape from the barrier, 
as can be observed from the vanishing escape rate at high 
friction. 


lim rjj/2 -h Cj — 2 




r. 


(31) 


In order for the largest positive unbound solution (the 
highest root of the characteristic equation whose solution 
gives Eq. (30)) to be the physical barrier frequency, the 


friction would have to be negative] Tjj = —uijfujb. 

For the (artificial) case of a reaction whose reaction co¬ 
ordinate is solely a single non-centroid normal mo de, the 
above result is corroborated by Kramers theorj)23MIin! 
which states that the transmission coefficient K{t) de¬ 
creases with friction as 


lim nit) ~ \/1 -\- — a 


(32) 


harmonic oscillator (where the external frequency must 
also be considered), and can be determined without knowl¬ 
edge of the frequencies present in the external potential. 
Obviously, chemical systems will not be purely harmonic 
but in many systems (such as vibrating bond) this will 
be a reasonable approximation. 

Previous literature has explored a range of scaled fric¬ 
tion matrices of AF and found A = 1/2 to be optimal for 
some spectra, justifying this on the grounds of optimal 
sampling of the harmonic ring polymer potential energjl^, 
but also finding there to be a wide range of A close to 
A = 1/2 in which results are broadly similar (as also 
seen in Fig. [^. We suspect that a numerically favourable 
value of A = 1/2 is due to interplay between shifting the 
frequencies of the higher normal modes to the external 
frequency (implying A = 1), moving the maximum in the 
spectral peak (implying A = 2“^^^), and avoiding harsh 
damping which would decorrelate the modes too quickly 
to capture their dynamics and broaden spectral peaka^Hl 
(implying the weakest possible friction which removes 
spurious resonances). We can certainly find no reason to 
use A > 1. 


where a = Tjj I2u]j and Cvj is the barrier frequency in 
ring-polymer space defined above. However, the trans¬ 
mission coefficient across a parabolic barrier is unity in 
Matsubara dynamics, and adding friction in TRPMD will 
only decrease this. Consequently, application of friction 
to RPMD will not ameliorate the qualitative problems 
with the RPMD higher normal modes at a barrier, and 
in some cases will worsen them. 


V. DISCUSSION 

The friction matrix F^y/ = 2\Cjj\5jji obtained sec¬ 
tion |IV A| corresponds to critical damping of the ring 
polymer springs in the absence of an external potential, 
but not critical damping of the ring polymer modes in a 


This definition of the friction matrix means that Too = 0 
(for all A) meaning that the centroid is unthermostatted, 
so all the results which have previously been derived 
for TRPMD, such as its short-time error compared to 
the quantum resullP^, still hold. A new result is that 
TRPMD, like RPMD, will have the exact Matsubara 
force on the centroid, since the error term does not act 
upon the centroid. Like RPMD and CMD but not LSC- 
IVR, the TRPMD dynamics will also satisfy detailed 
balanc^^, and the error scaling in the higher normal 
modes in time will be the same as that for RPMD, namely 
oc 1/ph. 

This choice of friction matrix also means that the 
TRPMD correlation function of a linear operator will 
deviate from the Matsubara correlation due to higher- 
order coupling between the centroid dynamics and the 
































damping (and random kicks) of the higher normal modes 
via anharmonicity in the potential. This causes slight 
broadening of spectral lines (a far smaller issue than the 
curvature problem of CMD or the spurious resonances of 
RPMD)P^, but the extra friction noticeably slows reaction 
rates beneath the crossover temperature^ where the un¬ 
bound and thermostatted higher normal modes are part 
of the optimal dividing surfac^^. For nonlinear operators, 
TRPMD (like RPMD) would be expected to break down 
faster than for linear operators due to the error term only 
acting directly on the higher normal modes, though the 
example of the position-squared autocorrelation function 
given above suggests that with a careful choice of friction 
the breakdown may not be too drastic. 

Although the analysis for a parabolic barrier in sec¬ 
tion |IV C| does not suggest that the TRPMD rate will 
ever be closer to the Matsubara (and therefore quan¬ 
tum) rate, TRPMD could be computationally advisable 
above the crossover temperature (where passage over the 
barrier is dominated by motion of the unthermostatted 
centroid) since the TRPMD trajectories may sample the 
path-integral phase space more efficiently than the RPMD 
trajectories^, and same may be true for other observ¬ 
able properties whic h are dominated by barrier crossing, 
such as diffusiorPSED xhe TRPMD time-evolution is also 
simpler computationally since the same dynamics can be 
used for thermostatting and computation of the correla¬ 
tion functiorPSl. Nevertheless, these results suggest that 
beneath the crossover temperature, TRPMD is not to be 
advised for reaction rates, a result broadly supported by 
numerical tests in one-dimensional and multidimensional 
gas-phase system^. 

All the results presented here generalize immediately 
to multidimensional systems, where the friction is applied 
in F{N — 1) normal modes and springs exist between 
N replicas of the physical system. For nonlinear opera¬ 
tors one cannot, in general, easily relate the Kubo and 
Generalized Kubo forms [Eqs. 0 and (|^] (the position- 
squared operator explored above being an exception). For 
reaction rates involving the highly nonlinear flux and side 
operators this is resolved by relating the generalized Kubo 
form to the exact quantum expression when there is no 
recrossing of the path-integral dividing surface (and those 
orthogonal to it in path-integral space) by the exact quan¬ 
tum dynamics of the systeirPlHlII^ such that TRPMD rate 
theory will give the exact quantum rate in the absence 
of recrossing by either the TRPMD dynamics or exact 
quantum dynamics.^ 


VI. CONCLUSIONS 

In this article we have shown, for the first time, how to 
obtain thermostatted ring polymer molecular dynamics 
(TRPMD) from exact quantum dynamics by a series of ap¬ 
proximations, each with an analytic error term. We firstly 
discard fluctuations of the highest N — M normal modes 
from the exact quantum time evolution, giving Matsub¬ 


ara dynamicJ^. To derive a computationally tractable 
approximation to Matsubara dynamics, we perform a con¬ 
tour integral in the momenta (where we assume the edge 
terms to be zero), giving a correlation function with the 
(real) ring polymer distribution, but whose Liouvillian is 
complex. We then replace the imaginary part of the com¬ 
plex Liouvillian with a white-noise Fokker-Planck term, 
giving TRPMD. 

Each of these approximations has its limitations and 
benefits. The primary consequence of discarding the 
fluctuations of the highest normal modes from the exact 
quantum dynamics (leading to Matsubara dynamics) is 
neglect of interference effects and mixing of quantum 
states. In physical systems this is seen as the failure of 
TRPMD to replicate the Fermi resonance in the Zundel 
catiorPEl (CMD and RPMD also fail her^^, as would be 
expected as they too are approximations to Matsubara 
dynamic^^. However, discarding these fluctuations leads 
to a classical-like dynamics which preserves the quantum 
Boltzmann distribution.^ 

We then show that a careful choice of the friction matrix 
(which is system independent and known in advance) will 
cause all ring polymer normal modes to oscillate at the 
correct frequency in a harmonic potential, and therefore 
will reproduce the correct frequency of oscillation of the 
position-squared autocorrelation function and the correct 
t = 0 value; neither CMD nor RPMD will replicate both 
of these properties. However, the oscillations’ amplitude 
is too s mall, and we suspect that a generalized Langevin 
equatiorf^J^ may be more successful than a simple white 
noise thermostat, in that it may be constructed to pro¬ 
duce the correct frequency of oscillation of the higher 
normal modes but with smaller damping (and maybe 
even the correct maxima in the position and momentum 
autocorrelation spectra)!^. The same analysis, but for 
a parabolic barrier, shows that no physical friction pa¬ 
rameter will solve the qualitative innacurracies in the 
higher ring polymer norma l mod e fluctuations. Usage of 
unphysical negative frictioip21£2l 3,3 a possible solution to 
this problem is left as further work. Future research could 
also include extension to non-adia batic systems where 

RPMD has been successful .1221541641161 

In closing, the results presented here give an a pri¬ 
ori prescription for when to use TRPMD: it should be 
used for computation of spectra and other properties of 
bound systems where the correct oscillation frequencies 
are required, and avoided for rate calculation beneath the 
crossover temperature. 


Appendix A: Matsubara modes 

The ring-polymer normal modes are defined as 

T.. 
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where j = —N/2 + 1,..., 0,..., N/2 and likewise for P, 
where 

r j = 0 

_ I ^/Wsin(27rzj7Af) 1 < j < iV/2 - 1 

u \ j = N/2 

[ y^cos(27rij7fV) -N/2 + 1 < j <-1 

(A2) 

where the j = N/2 mode is omitted if N is oddf^ The 
transformation is not unitary, but defined such that the 
normal modes converge in the iV —>■ oo limit. This leads 
to frequencies in the complex Boltzmann distribution of 


UJn = 


2 sin(j7r/7V) 
13 Nh 


(A3) 


which, for large N and finite j, become the Matsubara 
frequencies^ 


27rj 


a),- = lim oja = 

N^oo I3H 


(A4) 


The observables A(Q) and B(Q) are obtained by making 
by substituting 


(M-l)/2 

qi = ^ TijVNQj 

j=-{M-l)/2 


(AS) 


into A(q) and B{q) respectively, which also leads to a 
‘Matsubara potential’, 

N-l / (M-l)/2 \ 

= E (A6) 

i=0 J 


Appendix B: Equivalence of quantum and Matsubara 
correlation functions 


To show that the Matsubara correlation function is 


equivalent to Eq. (23), we firstly calculate the Matsubara 


correlation function using the harmonic analysis in the 
supplementary material of Ref. [371 giving 


cjSUt) = 




(M-l)/2 

cos{2uJht) E 


(M-l)/2 


E 


1 - {ojj/uJhY 


1 + (ujNuJh)'^ ^ (1 + 

j = -(M-l)/2 ^ i=-(M-l)/2 ^ \ 3 / H)) 


(M-l)/2 

E 


(M-l)/2 

E 


1 


1 


j = -(M-l)/2 fc=-(M-l)/2 


1 + {Cjj/uJh)'^ 1 + {uk/uJhY 


(Bl) 


The Matsubara frequency summatiorf^ is performed by 
examining the integral 


dz 


cot(z) 


(B2) 


around a circle of infinite radius, origin zero, we find 


(M-l)/2 

;coth(a:) = lim > 

' ' Z. -/ 


M-^-oo ' l + ( 77 r/x)^ 


(B3) 


and by differentiation of Eq. (B3), that 


Appendix C: Analyticity in the complex plane 

Consider an observable B(P, Q, t), which is propagated 
by the Liouvillian 


£ = (VpEf) • Vq - (Vqi?) • Vp 


(Cl) 


where H is the Hamiltonian of the system and an ana¬ 
lytic, but not necessarily real, function of P and Q. The 
propagation is formally 


-i?(P,Q,t) =£H(P,Q,t) 


i?(P,Q,t) =e^‘i?(P,Q,0) 


(C2) 

(C3) 


(M-l)/2 / ■ / N2 

7 coth^(x) - 1 ] = lim E 71- / A\2 ' 


(B4) 


This (obviously) requires B(P,Q,t) to be s ingle valued, 
and the exponentiated expression Eq. (C3| to exist. If 


B(P, Q, t) is an analytic function for all values of z, then 
(by the Cauchy-Riemann relations) 


Subsituting x = Phuih/2 into Eq. (B3) and Eq. (B4), and 


these expressions into Eq. (Bl) gives Eq. (23) as required. 


d 

W* 


H(P,Q,t)=0Vi 


(C4) 
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where P* is the complex conjugate of Pj (and likewise for 
Q*j)- If i? is analytic then 

= (C5) 


which means that (using H being continuous, Schwarz’ 
theorem and therefore the commuta¬ 

tion relations exist 


d 

W* 


L=L 


d 

W; 


(C6) 


Using the definition of an exponential as its power expan¬ 
sion we then see, 


^B(P,Q,t)=^e^‘B(P,Q,0) 

U-Tj 

= 0 


(C7) 

(C8) 

(C9) 


so B{P, Q,t) remains an analytic function of Pj for all 
time (and likewise for Qj). This is true Vj (and Vt), 
and by Hartog’s Theorem, true for i?(P, Q, t) everywhere. 
This means that B{P,Q,t) obeys the Cauchy Riemann 
relations and can have no poles in the complex plane. 
The Boltzmann distribution is also holomorphic, and 
provided that the zero-time observable A(P,Q,0) is also 
holomorphic (which almost all physical observables are) 
the entire integrand of Eq. Q will be. 

We then complete the square in the complex Matsub- 
ara distribution, giving Eq. where the edges of the 
rectangle used in the contour integration are 


f 

(M-l)/2 

pmujjQ — j 

/ dQ 

n ^ 

f dn, 

J 

II 

1 

1 

c 




lim 

TT—>■ —OO 


f 

(M-l)/2 

pmojjQ — j 

/ dQ 

n • 

f dn, 

J 

_i=-(M-l)/2 


(CIO) 


where iTj = ^Pj, Ilj = SP,-, and Matsubara 

Liouvillian Eq. Q continued into the complex plane. 


The edge terms can be proven to be zero in a number 
of limits. Specifically, for A(Q) and i?(Q) which are at 
most exponential in P and/or Q, the edge terms will 
vanish when the trajectories are real (II = 0) where 
conservation of energy arguments can be used in a bound 
system and in a scattering system whose potential tends 
to a constant value far out. The edges will also be zero in 
any system at t = 0 where the momentum integral can be 
evaluated analytically, and where discarding (and 
thereby keeping the trajectories real) is no approximation, 
namely up to 0{t^) f or nonlinear operators and 0(t®) for 

linear operator^^HSH 


For systems where the trajectories are known analyti¬ 
cally, such as a free particle, parabolic well and barrier, 
even though Tv{t) —)■ oo as 7r(0) —)■ oo, careful considera¬ 
tion of the limits and application of rHopital’s rule shows 
that the edge term still vanishes. 


Despite the above promising results, trajectories in 
the complex plane are frequently not bounded^ and in 
general it is difficult to determine whether or not terms 
of the form in Eq. (CIOI will converg^^for any general 
potential. A proof of whether S{t) can be neglected in 
any general case is left as further work. 


Appendix D: TRPMD position-squared correlation functions 


The correlation function can be evaluated by consider¬ 
ing each normal mode separately and derivin g the c orre- 
lation function for a single harmonic oscillatop32lIIZl. For 
0 < A < 1, the correlation function is 


''[^1 __ L 




(M-l)/2 

E 

3 = -{M-1)/2 


■’h +^j 


1 2e 

'To 

cos{LUjt) + sin(wjt) 

2 (M-l)/2 

1 V ’ 

1 Uji + U!^ 

^ cc/ + 1 





(Dl) 


with w, defined in Eq. (21). If A > 1, we define 


Jcut — 


phiuh 


27rVA2 - 1, 


(D2) 
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where [-J is the floor function, so all modes with |j| > jcut will be overdamped. The correlation function is 


(f) — ^ 

q^q-^y ! « 2^2 


E 


2 


J = -Jcut “ 

M/2-1 

E 


—r • f 

o 1 -J 


ujI + ujj 


1 + 


4C," 


rT 

_ 33 


4^2 I cos(2Cjt) 


sin(2^jt) 


(M-l)/2 
fc=-(M-l)/2 " 


1 


W? + W? 


—r t 


J=Jcut + l 




• 



2 I cosh(2(jt) + —^ sinh(2^jt) 

Ci 


(M-l)/2 

E 

/c=-(M-l)/2 


Ujt + Wt 


(D3) 


where we have noted that contributions from modes j 
and —j are the same. If a mode is critically damped 
then the term in square brackets for that mode becomes 

2 + r|//2 + 2r,4t. 
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